Nejjednodušší evoluční PDR (?) -- convekce

Transport veličiny $u$ podél toku o rychlosti $v$ je popsán rovnicí $$\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} = 0$$



In [40]:

    
%matplotlib inline
from pylab import *
from numpy import *



In [41]:

    
#dimensions of the computational domain:
maxx = 10.
maxt = 20.

#flow velocity
v = 0.5

We discretize the simulation domain: $$t_n = t_0 + n\Delta t$$ $$x_i = x_0 + i\Delta x$$ and seek for the numerical approximation $U^n_i$ to the actual solution $U(t_n, x_i)$ on the grid points



In [42]:

    
#discretization parameters
nx = 100   # number of unknown grid points in spatial direction
CFL = 1.0  # Courant constant v*dt/dx

def wave_init(maxx, maxt, v, nx, CFL):
    #choose time step according to CFL condition
    dx = maxx/(nx+1)
    dt = CFL*dx/v
    nt = int(maxt/dt)+1
    
    #define array for storing the solution
    U = zeros((nt, nx+2))
    
    x = arange(nx+2)*dx
    t = arange(nt)*dt
    return U, dx, dt, x, t

To obtain an explicit scheme for propagating the solution in time, we may replace the derivatives with forward differences in time and central differences in space (FTCS) $$\frac{U^{n+1}_i-U^n_i}{\Delta t} + c\frac{U^n_{i+1}-U^n_{i-1}}{2\Delta x}=0.$$ This leads to a simple explicit formula for $U^{n+1}_i$ $$U^{n+1}_i = U^n_i - \frac{c\Delta t}{2\Delta x}(U^n_{i+1}-U^n_{i-1})$$



In [43]:

    
def convection_solve(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,nt-1):
        U[it+1,xint] = U[it,xint] - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint-1])

Let's try to propagate a square initial condition



In [44]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
#U[0,:] = sin(x*5)
plot(U[0,:])
ylim(0,1.1)









    Out[44]:





(0, 1.1)



In [45]:

    
convection_solve(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[45]:





<matplotlib.collections.QuadMesh at 0x7f9ae56df978>

:-( Does it help if we choose smaller time step?



In [46]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[46]:





<matplotlib.collections.QuadMesh at 0x7f9ae52d1b00>

Not really, this solver fails to reproduce the expected solution. See Von Neumann stability analysis or matrix method...

The Lax(-Friedrichs) method stabilizes the solution by replacing $U^n_i$ with centered approximation $\frac{1}{2}(U^n_{i-1}+U^n_{i+1})$: $$U^{n+1}_i = \frac{1}{2}(U^n_{i-1}+U^n_{i+1}) - \frac{c\Delta t}{2\Delta x}(U^n_{i+1}-U^n_{i-1})$$



In [47]:

    
def convection_solve_Lax(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,nt-1):
        U[it+1,xint] = 0.5*(U[it,xint+1]+U[it,xint-1]) - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint-1])



In [49]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_Lax(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[49]:





<matplotlib.collections.QuadMesh at 0x7f9ae5219668>

The Courant-Friedrichs-Lewy (CFL) condition $$\frac{v\Delta t}{\Delta x}\leq 1$$ determines the stability. Look at the sensitivity to violating this condition:



In [51]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=1.05)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_Lax(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[51]:





<matplotlib.collections.QuadMesh at 0x7f9ae53e4e48>

What if we need better time resolution?



In [52]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_Lax(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[52]:





<matplotlib.collections.QuadMesh at 0x7f9ae5158860>

The stabilization also causes "damping" of the solution at smaller time steps.



In [53]:

    
def convection_solve_upwind(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,nt-1):
        U[it+1,xint] = U[it,xint] - v*dt/dx*(U[it,xint] - U[it, xint-1])



In [59]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_upwind(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[59]:





<matplotlib.collections.QuadMesh at 0x7f9ae1fc9ac8>



In [55]:

    
def convection_solve_LeapFrog(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    #start the solution with upwind method
    U[1,xint] = U[0,xint] - v*dt/dx*(U[0,xint] - U[0, xint-1])
    for it in range(1, nt-1):
        U[it+1,xint] = U[it-1,xint] - v*dt/dx*(U[it,xint+1] - U[it, xint-1])

U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LeapFrog(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[55]:





<matplotlib.collections.QuadMesh at 0x7f9ae509bcc0>



In [60]:

    
def convection_solve_LeapFrog(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    #start the solution with upwind method
    U[1,xint] = U[0,xint] - v*dt/dx*(U[0,xint] - U[0, xint-1])
    for it in range(1, nt-1):
        U[it+1,xint] = U[it-1,xint] - v*dt/dx*(U[it,xint+1] - U[it, xint-1])
                
        #extrapolate values on the outflow boundary
        U[it+1,-1] = U[it+1, -2]
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LeapFrog(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[60]:





<matplotlib.collections.QuadMesh at 0x7f9ae201d358>



In [61]:

    
def convection_solve_LeapFrog(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    #start and smooth the solution with upwind method
    for it in range(0,min(50, nt-1)):
        U[it+1,xint] = U[it,xint] - v*dt/dx*(U[it,xint] - U[it, xint-1])
        #U[it+1,xint] = 0.5*(U[it,xint+1]+U[it,xint-1]) - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint-1])
    for it in range(min(50, nt-1),nt-1):
        U[it+1,xint] = U[it-1,xint] - v*dt/dx*(U[it,xint+1] - U[it, xint-1])
                
        #extrapolate values on the outflow boundary
        U[it+1,-1] = U[it+1, -2]*2 - U[it+1, -3]
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LeapFrog(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)









    Out[61]:





<matplotlib.collections.QuadMesh at 0x7f9ae1fcf908>



In [ ]:

    
def convection_solve_LaxWendroff(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,nt-1):
        # less efficient but more readable:
        Uplus  = 0.5*(U[it,xint+1] + U[it, xint]) - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint])
        Uminus = 0.5*(U[it,xint] + U[it, xint-1]) - 0.5*v*dt/dx*(U[it,xint] - U[it, xint-1])
        U[it+1,xint] = U[it,xint] - v*dt/dx*(Uplus - Uminus)



In [ ]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LaxWendroff(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)



In [ ]:

    
def convection_solve_LaxWendroff(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,nt-1):
        # less efficient but more readable:
        Uplus  = 0.5*(U[it,xint+1] + U[it, xint]) - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint])
        Uminus = 0.5*(U[it,xint] + U[it, xint-1]) - 0.5*v*dt/dx*(U[it,xint] - U[it, xint-1])
        U[it+1,xint] = U[it,xint] - v*dt/dx*(Uplus - Uminus)
        
        #extrapolate values on the outflow boundary
        U[it+1,-1] = U[it+1, -2]*2 - U[it+1,-3]



In [ ]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LaxWendroff(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)



In [ ]:

    
def convection_solve_LaxWendroff(U, dx, dt, nx, nt):
    xint = arange(1, nx+1)
    for it in range(0,min(50, nt-1)):
        U[it+1,xint] = U[it,xint] - v*dt/dx*(U[it,xint] - U[it, xint-1])
    for it in range(min(50, nt-1),nt-1):
        # less efficient but more readable:
        Uplus  = 0.5*(U[it,xint+1] + U[it, xint]) - 0.5*v*dt/dx*(U[it,xint+1] - U[it, xint])
        Uminus = 0.5*(U[it,xint] + U[it, xint-1]) - 0.5*v*dt/dx*(U[it,xint] - U[it, xint-1])
        U[it+1,xint] = U[it,xint] - v*dt/dx*(Uplus - Uminus)
        
        #extrapolate values on the outflow boundary
        U[it+1,-1] = U[it+1, -2]*2 - U[it+1,-3]

U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=.1)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_LaxWendroff(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)

Can we break the CFL condition with some naive implementation of implicit scheme?



In [ ]:

    
from scipy.sparse import dia_matrix, eye
from scipy.sparse.linalg import factorized
def d1matrix(nelem):
    elements = ones((3,nelem))
    elements[1,:] *= 0
    elements[0,:] *= -1
    return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()

def convection_solve_CN(U, dx, dt, nx, nt):
    alpha = -v*dt/(dx*4.)
    M1 = eye(nx)-d1matrix(nx)*alpha
    M2 = eye(nx)+d1matrix(nx)*alpha
    LU = factorized(M1.tocsc())
    for it in range(0,nt-1):
        U[it+1,1:-1] = LU(M2.dot(U[it,1:-1]))
        
        #extrapolate values on the outflow boundary
        U[it+1,-1] = U[it+1, -2]



In [ ]:

    
U, dx, dt, x, t = wave_init(maxx, maxt, v, nx, CFL=2)
U[0,:] = 0.
U[0, (x<maxx*0.4) & (x>maxx*0.2)] = 1.
convection_solve_CN(U, dx, dt, nx, len(t))
pcolormesh(U, rasterized=True, vmin=-2, vmax=2)

It sort of works, but not very well...



In [ ]: